# Comprehensive SEM Analysis with Multiple Competing Models
# Testing original hypotheses plus alternative and reversed causal models

# ================= SECTION 1: SETUP AND DATA PREPARATION =================
# Load necessary packages
library(lavaan)     # For SEM analysis
library(psych)      # For descriptive statistics and reliability analysis
library(semTools)   # For additional SEM tools
library(dplyr)      # For data manipulation

# Import data
data <- read.csv("Charcter data.csv", header = TRUE, stringsAsFactors = FALSE)

# Convert all relevant variables to numeric
for(col in names(data)) {
  if(!is.numeric(data[[col]])) {
    cat("Converting", col, "to numeric\n")
    data[[col]] <- as.numeric(as.character(data[[col]]))
  }
}

# Check for missing values
missing_summary <- colSums(is.na(data))
cat("Variables with missing values:\n")
print(missing_summary[missing_summary > 0])

# ================= SECTION 2: IMPROVED MEASUREMENT MODELS =================

# Character Strength measurement model - this had good fit, so keep as is
char_strength_model <- '
  # Character Strength as a latent variable with 6 virtue indicators
  CharacterStrength =~ WISDOM + COURAGE + HUMANITY + JUSTICE + TEMPERANCE + TRANSCENDENCE
'

# PROPERLY CORRECTED Well-being model - using a simpler approach
wellbeing_model <- '
  # Well-being as a latent variable with 5 PERMA indicators
  Wellbeing =~ Positive_Negative_emotions + Engagement + Relationships + Meaning + Accomplishment
  
  # Adding more flexibility with additional correlated errors
  Positive_Negative_emotions ~~ Engagement
  Positive_Negative_emotions ~~ Relationships
  Engagement ~~ Accomplishment
  Relationships ~~ Meaning
  Meaning ~~ Accomplishment
'

# PROPERLY CORRECTED Affective Temperament model 
# This model had serious problems, so we'll use a more conservative approach
temperament_model <- '
  # Affective Temperament as a latent variable with 4 indicators
  # Removed Hyperthymic (was non-significant) and will use fewer correlated errors
  AffectiveTemperament =~ Cyclothymic + Depressive + Irritable + Anxious
  
  # Allow only one correlated error instead of multiple
  Cyclothymic ~~ Irritable
'

# ================= SECTION 3: FITTING MEASUREMENT MODELS =================

# Fit Character Strength measurement model
tryCatch({
  fit_char_strength <- cfa(char_strength_model, data = data, missing = "fiml")
  summary_char_strength <- summary(fit_char_strength, fit.measures = TRUE, standardized = TRUE)
  cat("\nCHARACTER STRENGTH MEASUREMENT MODEL:\n")
  print(summary_char_strength)
}, error = function(e) {
  cat("Error fitting Character Strength model:", e$message, "\n")
})

# Fit improved Well-being measurement model
tryCatch({
  fit_wellbeing <- cfa(wellbeing_model, data = data, missing = "fiml")
  summary_wellbeing <- summary(fit_wellbeing, fit.measures = TRUE, standardized = TRUE)
  cat("\nWELL-BEING MEASUREMENT MODEL (IMPROVED):\n")
  print(summary_wellbeing)
}, error = function(e) {
  cat("Error fitting Well-being model:", e$message, "\n")
})

# Fit improved Temperament measurement model
tryCatch({
  fit_temperament <- cfa(temperament_model, data = data, missing = "fiml")
  summary_temperament <- summary(fit_temperament, fit.measures = TRUE, standardized = TRUE)
  cat("\nAFFECTIVE TEMPERAMENT MEASUREMENT MODEL (IMPROVED):\n")
  print(summary_temperament)
}, error = function(e) {
  cat("Error fitting Affective Temperament model:", e$message, "\n")
})

# ================= SECTION 4: STRUCTURAL MODELS (ORIGINAL HYPOTHESES) =================

# Direct Effects Model (H1) - with improved measurement models
direct_model <- '
  # Measurement models
  CharacterStrength =~ WISDOM + COURAGE + HUMANITY + JUSTICE + TEMPERANCE + TRANSCENDENCE
  
  # Improved Well-being model
  Wellbeing =~ Positive_Negative_emotions + Engagement + Relationships + Meaning + Accomplishment
  Positive_Negative_emotions ~~ Engagement
  Positive_Negative_emotions ~~ Relationships
  Engagement ~~ Accomplishment
  Relationships ~~ Meaning
  Meaning ~~ Accomplishment
  
  # Structural model - direct path
  Wellbeing ~ c*CharacterStrength
'

# Full Mediation Model (H3) - with improved measurement models
full_mediation_model <- '
  # Measurement models
  CharacterStrength =~ WISDOM + COURAGE + HUMANITY + JUSTICE + TEMPERANCE + TRANSCENDENCE
  
  # Improved Affective Temperament model
  AffectiveTemperament =~ Cyclothymic + Depressive + Irritable + Anxious
  Cyclothymic ~~ Irritable
  
  # Improved Well-being model
  Wellbeing =~ Positive_Negative_emotions + Engagement + Relationships + Meaning + Accomplishment
  Positive_Negative_emotions ~~ Engagement
  Positive_Negative_emotions ~~ Relationships
  Engagement ~~ Accomplishment
  Relationships ~~ Meaning
  Meaning ~~ Accomplishment
  
  # Structural model - full mediation (no direct path)
  AffectiveTemperament ~ a*CharacterStrength
  Wellbeing ~ b*AffectiveTemperament
  
  # Define indirect effect
  indirect := a*b
'

# Partial Mediation Model (H2) - with improved measurement models
partial_mediation_model <- '
  # Measurement models
  CharacterStrength =~ WISDOM + COURAGE + HUMANITY + JUSTICE + TEMPERANCE + TRANSCENDENCE
  
  # Improved Affective Temperament model
  AffectiveTemperament =~ Cyclothymic + Depressive + Irritable + Anxious
  Cyclothymic ~~ Irritable
  
  # Improved Well-being model
  Wellbeing =~ Positive_Negative_emotions + Engagement + Relationships + Meaning + Accomplishment
  Positive_Negative_emotions ~~ Engagement
  Positive_Negative_emotions ~~ Relationships
  Engagement ~~ Accomplishment
  Relationships ~~ Meaning
  Meaning ~~ Accomplishment
  
  # Structural model - partial mediation
  AffectiveTemperament ~ a*CharacterStrength
  Wellbeing ~ c_prime*CharacterStrength + b*AffectiveTemperament
  
  # Define direct, indirect, and total effects
  indirect := a*b
  total := c_prime + (a*b)
'

# ================= SECTION 5: ALTERNATIVE THEORETICAL MODELS =================

# Reversed Causal Model - tests opposite direction of causality
reversed_causal_model <- '
  # Measurement models
  CharacterStrength =~ WISDOM + COURAGE + HUMANITY + JUSTICE + TEMPERANCE + TRANSCENDENCE
  
  # Improved Affective Temperament model
  AffectiveTemperament =~ Cyclothymic + Depressive + Irritable + Anxious
  Cyclothymic ~~ Irritable
  
  # Improved Well-being model
  Wellbeing =~ Positive_Negative_emotions + Engagement + Relationships + Meaning + Accomplishment
  Positive_Negative_emotions ~~ Engagement
  Positive_Negative_emotions ~~ Relationships
  Engagement ~~ Accomplishment
  Relationships ~~ Meaning
  Meaning ~~ Accomplishment
  
  # Structural model - reversed causality
  AffectiveTemperament ~ a*Wellbeing
  CharacterStrength ~ c_prime*Wellbeing + b*AffectiveTemperament
  
  # Define direct, indirect, and total effects
  indirect := a*b
  total := c_prime + (a*b)
'

# Alternative Model - Affective Temperament as antecedent
alternative_model <- '
  # Measurement models
  CharacterStrength =~ WISDOM + COURAGE + HUMANITY + JUSTICE + TEMPERANCE + TRANSCENDENCE
  
  # Improved Affective Temperament model
  AffectiveTemperament =~ Cyclothymic + Depressive + Irritable + Anxious
  Cyclothymic ~~ Irritable
  
  # Improved Well-being model
  Wellbeing =~ Positive_Negative_emotions + Engagement + Relationships + Meaning + Accomplishment
  Positive_Negative_emotions ~~ Engagement
  Positive_Negative_emotions ~~ Relationships
  Engagement ~~ Accomplishment
  Relationships ~~ Meaning
  Meaning ~~ Accomplishment
  
  # Structural model - alternative causal sequence
  CharacterStrength ~ a*AffectiveTemperament
  Wellbeing ~ c_prime*AffectiveTemperament + b*CharacterStrength
  
  # Define direct, indirect, and total effects
  indirect := a*b
  total := c_prime + (a*b)
'

# ================= SECTION 6: FITTING ALL STRUCTURAL MODELS =================

# Fit Direct Effects Model
tryCatch({
  fit_direct <- sem(direct_model, data = data, missing = "fiml")
  summary_direct <- summary(fit_direct, fit.measures = TRUE, standardized = TRUE)
  cat("\nDIRECT EFFECTS MODEL:\n")
  print(summary_direct)
}, error = function(e) {
  cat("Error fitting Direct Effects model:", e$message, "\n")
})

# Fit Full Mediation Model
tryCatch({
  fit_full_mediation <- sem(full_mediation_model, data = data, missing = "fiml")
  summary_full_mediation <- summary(fit_full_mediation, fit.measures = TRUE, standardized = TRUE)
  cat("\nFULL MEDIATION MODEL:\n")
  print(summary_full_mediation)
}, error = function(e) {
  cat("Error fitting Full Mediation model:", e$message, "\n")
})

# Fit Partial Mediation Model
tryCatch({
  fit_partial_mediation <- sem(partial_mediation_model, data = data, missing = "fiml")
  summary_partial_mediation <- summary(fit_partial_mediation, fit.measures = TRUE, standardized = TRUE)
  cat("\nPARTIAL MEDIATION MODEL:\n")
  print(summary_partial_mediation)
}, error = function(e) {
  cat("Error fitting Partial Mediation model:", e$message, "\n")
})

# Fit Reversed Causal Model
tryCatch({
  fit_reversed_causal <- sem(reversed_causal_model, data = data, missing = "fiml")
  summary_reversed_causal <- summary(fit_reversed_causal, fit.measures = TRUE, standardized = TRUE)
  cat("\nREVERSED CAUSAL MODEL (Well-being → Affective Temperament → Character Strengths):\n")
  print(summary_reversed_causal)
}, error = function(e) {
  cat("Error fitting Reversed Causal model:", e$message, "\n")
})

# Fit Alternative Model
tryCatch({
  fit_alternative <- sem(alternative_model, data = data, missing = "fiml")
  summary_alternative <- summary(fit_alternative, fit.measures = TRUE, standardized = TRUE)
  cat("\nALTERNATIVE MODEL (Affective Temperament → Character Strengths → Well-being):\n")
  print(summary_alternative)
}, error = function(e) {
  cat("Error fitting Alternative model:", e$message, "\n")
})

# ================= SECTION 7: COMPARING ORIGINAL MODELS =================

# Compare the original three models if they were successfully fitted
if(exists("fit_direct") && exists("fit_full_mediation") && exists("fit_partial_mediation")) {
  # Compare fit indices
  fit_measures <- data.frame(
    Model = c("Direct Effects", "Full Mediation", "Partial Mediation"),
    ChiSq = c(fitMeasures(fit_direct, "chisq"), 
              fitMeasures(fit_full_mediation, "chisq"),
              fitMeasures(fit_partial_mediation, "chisq")),
    df = c(fitMeasures(fit_direct, "df"),
           fitMeasures(fit_full_mediation, "df"),
           fitMeasures(fit_partial_mediation, "df")),
    CFI = c(fitMeasures(fit_direct, "cfi"),
            fitMeasures(fit_full_mediation, "cfi"),
            fitMeasures(fit_partial_mediation, "cfi")),
    TLI = c(fitMeasures(fit_direct, "tli"),
            fitMeasures(fit_full_mediation, "tli"),
            fitMeasures(fit_partial_mediation, "tli")),
    RMSEA = c(fitMeasures(fit_direct, "rmsea"),
              fitMeasures(fit_full_mediation, "rmsea"),
              fitMeasures(fit_partial_mediation, "rmsea")),
    SRMR = c(fitMeasures(fit_direct, "srmr"),
             fitMeasures(fit_full_mediation, "srmr"),
             fitMeasures(fit_partial_mediation, "srmr")),
    AIC = c(fitMeasures(fit_direct, "aic"),
            fitMeasures(fit_full_mediation, "aic"),
            fitMeasures(fit_partial_mediation, "aic")),
    BIC = c(fitMeasures(fit_direct, "bic"),
            fitMeasures(fit_full_mediation, "bic"),
            fitMeasures(fit_partial_mediation, "bic"))
  )
  
  cat("\nORIGINAL MODEL FIT COMPARISON:\n")
  print(fit_measures)
  
  # Compare Partial vs Full Mediation (nested models)
  anova_result <- anova(fit_full_mediation, fit_partial_mediation)
  cat("\nChi-Square Difference Test (Full vs. Partial Mediation):\n")
  print(anova_result)
  
  # Interpret the result
  if(anova_result$`Pr(>Chisq)`[2] < 0.05) {
    cat("\nThe chi-square difference test is significant (p <", anova_result$`Pr(>Chisq)`[2], ").")
    cat("\nThis suggests that the Partial Mediation model fits significantly better than the Full Mediation model.")
    cat("\nHypothesis H3 (full mediation) is NOT supported.\n")
  } else {
    cat("\nThe chi-square difference test is NOT significant (p =", anova_result$`Pr(>Chisq)`[2], ").")
    cat("\nThis suggests that the Full Mediation model is not significantly worse than the Partial Mediation model.")
    cat("\nAccording to the principle of parsimony, the Full Mediation model might be preferred.")
    cat("\nHypothesis H3 (full mediation) is supported.\n")
  }
}

# ================= SECTION 8: COMPREHENSIVE MODEL COMPARISON =================

# Compare all five models if they were successfully fitted
if(exists("fit_direct") && 
   exists("fit_full_mediation") && 
   exists("fit_partial_mediation") && 
   exists("fit_reversed_causal") && 
   exists("fit_alternative")) {
  
  # Compare fit indices
  fit_measures_all <- data.frame(
    Model = c("Direct Effects", "Full Mediation", "Partial Mediation", 
              "Reversed Causal", "Alternative"),
    ChiSq = c(fitMeasures(fit_direct, "chisq"), 
              fitMeasures(fit_full_mediation, "chisq"),
              fitMeasures(fit_partial_mediation, "chisq"),
              fitMeasures(fit_reversed_causal, "chisq"),
              fitMeasures(fit_alternative, "chisq")),
    df = c(fitMeasures(fit_direct, "df"),
           fitMeasures(fit_full_mediation, "df"),
           fitMeasures(fit_partial_mediation, "df"),
           fitMeasures(fit_reversed_causal, "df"),
           fitMeasures(fit_alternative, "df")),
    CFI = c(fitMeasures(fit_direct, "cfi"),
            fitMeasures(fit_full_mediation, "cfi"),
            fitMeasures(fit_partial_mediation, "cfi"),
            fitMeasures(fit_reversed_causal, "cfi"),
            fitMeasures(fit_alternative, "cfi")),
    TLI = c(fitMeasures(fit_direct, "tli"),
            fitMeasures(fit_full_mediation, "tli"),
            fitMeasures(fit_partial_mediation, "tli"),
            fitMeasures(fit_reversed_causal, "tli"),
            fitMeasures(fit_alternative, "tli")),
    RMSEA = c(fitMeasures(fit_direct, "rmsea"),
              fitMeasures(fit_full_mediation, "rmsea"),
              fitMeasures(fit_partial_mediation, "rmsea"),
              fitMeasures(fit_reversed_causal, "rmsea"),
              fitMeasures(fit_alternative, "rmsea")),
    SRMR = c(fitMeasures(fit_direct, "srmr"),
             fitMeasures(fit_full_mediation, "srmr"),
             fitMeasures(fit_partial_mediation, "srmr"),
             fitMeasures(fit_reversed_causal, "srmr"),
             fitMeasures(fit_alternative, "srmr")),
    AIC = c(fitMeasures(fit_direct, "aic"),
            fitMeasures(fit_full_mediation, "aic"),
            fitMeasures(fit_partial_mediation, "aic"),
            fitMeasures(fit_reversed_causal, "aic"),
            fitMeasures(fit_alternative, "aic")),
    BIC = c(fitMeasures(fit_direct, "bic"),
            fitMeasures(fit_full_mediation, "bic"),
            fitMeasures(fit_partial_mediation, "bic"),
            fitMeasures(fit_reversed_causal, "bic"),
            fitMeasures(fit_alternative, "bic"))
  )
  
  cat("\nCOMPREHENSIVE MODEL FIT COMPARISON:\n")
  print(fit_measures_all)
  
  # Identify best fitting model based on multiple criteria
  best_aic <- which.min(fit_measures_all$AIC)
  best_bic <- which.min(fit_measures_all$BIC)
  best_cfi <- which.max(fit_measures_all$CFI)
  best_rmsea <- which.min(fit_measures_all$RMSEA)
  best_srmr <- which.min(fit_measures_all$SRMR)
  
  cat("\nMODEL SELECTION SUMMARY:\n")
  cat("Best model by AIC: ", fit_measures_all$Model[best_aic], "\n")
  cat("Best model by BIC: ", fit_measures_all$Model[best_bic], "\n")
  cat("Best model by CFI: ", fit_measures_all$Model[best_cfi], "\n")
  cat("Best model by RMSEA: ", fit_measures_all$Model[best_rmsea], "\n")
  cat("Best model by SRMR: ", fit_measures_all$Model[best_srmr], "\n")
  
  # Compare alternative causal directions with direct AIC/BIC differences
  cat("\nALTERNATIVE MODEL COMPARISON (AIC & BIC DIFFERENCES):\n")
  
  aic_diff_rev <- fit_measures_all$AIC[3] - fit_measures_all$AIC[4]  # Partial - Reversed
  aic_diff_alt <- fit_measures_all$AIC[3] - fit_measures_all$AIC[5]  # Partial - Alternative
  
  bic_diff_rev <- fit_measures_all$BIC[3] - fit_measures_all$BIC[4]  # Partial - Reversed
  bic_diff_alt <- fit_measures_all$BIC[3] - fit_measures_all$BIC[5]  # Partial - Alternative
  
  cat("Partial Mediation vs. Reversed Causal AIC difference: ", aic_diff_rev, 
      "(negative means Partial Mediation is better)\n")
  cat("Partial Mediation vs. Alternative AIC difference: ", aic_diff_alt, 
      "(negative means Partial Mediation is better)\n")
  
  cat("Partial Mediation vs. Reversed Causal BIC difference: ", bic_diff_rev, 
      "(negative means Partial Mediation is better)\n")
  cat("Partial Mediation vs. Alternative BIC difference: ", bic_diff_alt, 
      "(negative means Partial Mediation is better)\n")
}

# ================= SECTION 9: DETAILED PARAMETER ANALYSIS =================

# Analyze the parameters of the best-fitting model (partial mediation)
if(exists("fit_partial_mediation")) {
  # Extract parameter estimates
  params <- parameterEstimates(fit_partial_mediation, standardized = TRUE)
  
  # Function to find parameters more reliably
  get_path_value <- function(params_df, lhs_val, rhs_val, op_val="~") {
    row_idx <- which(params_df$lhs == lhs_val & 
                       params_df$rhs == rhs_val & 
                       params_df$op == op_val)
    
    if(length(row_idx) > 0) {
      return(params_df$est.std[row_idx])
    } else {
      return(NA)
    }
  }
  
  # Extract path coefficients
  a_path <- get_path_value(params, "AffectiveTemperament", "CharacterStrength")
  b_path <- get_path_value(params, "Wellbeing", "AffectiveTemperament")
  c_prime <- get_path_value(params, "Wellbeing", "CharacterStrength")
  
  # Extract indirect and total effects
  indirect_effect <- params$est.std[params$lhs == "indirect" & params$op == ":="]
  total_effect <- params$est.std[params$lhs == "total" & params$op == ":="]
  
  # Display results
  cat("\nPARTIAL MEDIATION MODEL PARAMETER ESTIMATES:\n")
  cat("Path a (CharacterStrength → AffectiveTemperament):", round(a_path, 3), "\n")
  cat("Path b (AffectiveTemperament → Wellbeing):", round(b_path, 3), "\n")
  cat("Direct effect c' (CharacterStrength → Wellbeing):", round(c_prime, 3), "\n")
  cat("Indirect effect (a*b):", round(indirect_effect, 3), "\n")
  cat("Total effect (c' + a*b):", round(total_effect, 3), "\n")
  
  # Calculate proportion mediated
  if(!is.na(indirect_effect) && !is.na(total_effect) && total_effect != 0) {
    prop_mediated <- indirect_effect / total_effect
    cat("Proportion mediated:", round(prop_mediated, 3), "\n")
  }
}

# Do the same for the alternative model if it exists
if(exists("fit_alternative")) {
  # Extract parameter estimates
  params_alt <- parameterEstimates(fit_alternative, standardized = TRUE)
  
  # Extract path coefficients using the helper function
  a_path_alt <- get_path_value(params_alt, "CharacterStrength", "AffectiveTemperament")
  b_path_alt <- get_path_value(params_alt, "Wellbeing", "CharacterStrength")
  c_prime_alt <- get_path_value(params_alt, "Wellbeing", "AffectiveTemperament")
  
  # Extract indirect and total effects
  indirect_effect_alt <- params_alt$est.std[params_alt$lhs == "indirect" & params_alt$op == ":="]
  total_effect_alt <- params_alt$est.std[params_alt$lhs == "total" & params_alt$op == ":="]
  
  # Display results
  cat("\nALTERNATIVE MODEL PARAMETER ESTIMATES:\n")
  cat("Path a (AffectiveTemperament → CharacterStrength):", round(a_path_alt, 3), "\n")
  cat("Path b (CharacterStrength → Wellbeing):", round(b_path_alt, 3), "\n")
  cat("Direct effect c' (AffectiveTemperament → Wellbeing):", round(c_prime_alt, 3), "\n")
  cat("Indirect effect (a*b):", round(indirect_effect_alt, 3), "\n")
  cat("Total effect (c' + a*b):", round(total_effect_alt, 3), "\n")
  
  # Calculate proportion mediated
  if(!is.na(indirect_effect_alt) && !is.na(total_effect_alt) && total_effect_alt != 0) {
    prop_mediated_alt <- indirect_effect_alt / total_effect_alt
    cat("Proportion mediated:", round(prop_mediated_alt, 3), "\n")
  }
}

# Same for reversed causal model
if(exists("fit_reversed_causal")) {
  # Extract parameter estimates
  params_rev <- parameterEstimates(fit_reversed_causal, standardized = TRUE)
  
  # Extract path coefficients
  a_path_rev <- get_path_value(params_rev, "AffectiveTemperament", "Wellbeing")
  b_path_rev <- get_path_value(params_rev, "CharacterStrength", "AffectiveTemperament")
  c_prime_rev <- get_path_value(params_rev, "CharacterStrength", "Wellbeing")
  
  # Extract indirect and total effects
  indirect_effect_rev <- params_rev$est.std[params_rev$lhs == "indirect" & params_rev$op == ":="]
  total_effect_rev <- params_rev$est.std[params_rev$lhs == "total" & params_rev$op == ":="]
  
  # Display results
  cat("\nREVERSED CAUSAL MODEL PARAMETER ESTIMATES:\n")
  cat("Path a (Wellbeing → AffectiveTemperament):", round(a_path_rev, 3), "\n")
  cat("Path b (AffectiveTemperament → CharacterStrength):", round(b_path_rev, 3), "\n")
  cat("Direct effect c' (Wellbeing → CharacterStrength):", round(c_prime_rev, 3), "\n")
  cat("Indirect effect (a*b):", round(indirect_effect_rev, 3), "\n")
  cat("Total effect (c' + a*b):", round(total_effect_rev, 3), "\n")
  
  # Calculate proportion mediated
  if(!is.na(indirect_effect_rev) && !is.na(total_effect_rev) && total_effect_rev != 0) {
    prop_mediated_rev <- indirect_effect_rev / total_effect_rev
    cat("Proportion mediated:", round(prop_mediated_rev, 3), "\n")
  }
}

# ================= SECTION 10: BOOTSTRAP ANALYSIS (OPTIONAL) =================
# This is computationally intensive, so it's kept as an optional section

run_bootstrap <- FALSE  # Set to TRUE to run bootstrap analysis

if(run_bootstrap && exists("fit_partial_mediation")) {
  set.seed(123)  # For reproducibility
  
  cat("\nRunning bootstrap analysis for Partial Mediation model (this may take some time)...\n")
  
  tryCatch({
    # Bootstrap with 1000 samples
    boot_partial <- semTools::bootstrapLavaan(
      fit_partial_mediation, 
      R = 1000,
      type = "bca.simple",
      FUN = function(x) {
        est <- parameterEstimates(x, standardized = TRUE)
        
        # Extract parameters of interest
        a_path <- est$est.std[est$lhs == "AffectiveTemperament" & 
                                est$rhs == "CharacterStrength" & est$op == "~"]
        b_path <- est$est.std[est$lhs == "Wellbeing" & 
                                est$rhs == "AffectiveTemperament" & est$op == "~"]
        c_prime <- est$est.std[est$lhs == "Wellbeing" & 
                                 est$rhs == "CharacterStrength" & est$op == "~"]
        indirect <- est$est.std[est$lhs == "indirect" & est$op == ":="]
        total <- est$est.std[est$lhs == "total" & est$op == ":="]
        
        return(c(a_path = a_path, b_path = b_path, c_prime = c_prime, 
                 indirect = indirect, total = total))
      }
    )
    
    # Print bootstrap results
    cat("\nBOOTSTRAP RESULTS FOR PARTIAL MEDIATION MODEL:\n")
    print(boot_partial)
    
    # Calculate and display 95% confidence intervals
    boot_ci <- apply(boot_partial, 2, function(x) quantile(x, c(0.025, 0.975), na.rm = TRUE))
    cat("\n95% Bootstrap Confidence Intervals:\n")
    print(boot_ci)
    
  }, error = function(e) {
    cat("Error in bootstrap analysis:", e$message, "\n")
  })
}